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0^ ■ ABSTRACT 

^: 

Ph' Two dimensional realizations of self-consistent models for the "perfect elliptic 

O ' disks" were tested for global stability by gravitational N-body integration. 

' The family of perfect elliptic disk potentials have two isolating integrals; time 

■ independent distribution functions f{E, I2) which self-consistently reproduce 

■ the density distribution can be found numerically, using a modified marching 
^ ' scheme to compute the relative contributions of each member in a library 

c3 ■ of orbits. The possible solutions are not unique: for a given ellipticity, the 

models can have a range of angular momenta. Here results are presented for 
cases with minimal angular momentum, hence maximal random motion. As 
in previous work, N-body realizations were constructed using a modified quiet 
start technique to place particles on these orbits uniformly in action-angle space, 
making the initial conditions as smooth as possible. The most elliptical models 
initially showed bending instabilities; by the end of the run they had become 
slightly rounder. The most nearly axisymmetric models tended to become 
more elongated, reminiscent of the radial orbit instability in spherical systems. 
Between these extremes, there is a range of axial ratios 0.305 ^b/a ^ 0.570 for 
which the minimum streaming models appear to be stable. 
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1. Introduction 

Recent studies of elliptical galaxies and of bulges of spiral galaxies indicate that their 
figures are likely to be at least slightly triaxial (for reviews see Pinney 1982| ; |de Zeeuw fc 



Franx 1991 ; Bertin fc Stiavelli 1993|) . Most elliptical galaxies appear to be supported at 



least in part by anisotropics in the velocity distributions rather than by rapid rotation: see 
for example, the work on the dwarf elliptical galaxies NGC 147, 185 and 205 by [Bender 



Paquet fc Nieto (1991)| and [Held et al. (1992]|. A class of non-rotating potentials, known 



as the perfect ellipsoids, has been advanced as a possible model for elliptical galaxies (e.g. 
[de Zeeuw 1985[) . In these potentials, the mass density is stratified on concentric, similar 
ellipsoids, and is non-singular in the center. Many of the properties of these potentials can 
be derived analytically; the orbits all have three isolating integrals, and hence properties 
such as the time-averaged density distribution can be computed exactly. This simplifies 
the task of finding self- consistent models: time-steady phase-space distribution functions 
/(x, v) such that the resulting mass density generates the desired gravitational potential. 
[Statler (1987)[ and [leuben (1987)[ have demonstrated that distribution functions for the 



perfect ellipsoids, and the analogous two-dimensional elliptic disks, can be constructed. 
Various sub-families of the axisymmetric perfect ellipsoids have been tested for stability 
( [de Zeeuw fc Schwarzschild 19891 ; [Merritt fc Stiavelli 1991 [Merritt fc Hernquist 199ll ; 



[Robi.in fc de Zeeuw 1991] ). 

Flattened perfect ellipsoids could also be viewed as models for galactic bars. The only 
analytical bar models are Freeman's ( [19664 0! 0) bars, which are based upon a rotating 
two dimensional harmonic oscillator potential, and the perfect elliptic disk models, which 
have no figure rotation. [Iremaine fc de Zeeuw (1987)| showed that in the limit of the needle 



{b 0) the two dimensional perfect elliptic disk is neutrally stable. Prompted by the large 
streaming velocities seen in barred spiral galaxies, the stability of perfect elliptic disks 
with maximum angular momentum has already been studied ( P^jcvine fc Sparke 1994, LS[) . 
The roundest disks were unstable to spiral mode formation, the most elongated elliptical 
models were unstable to bending modes, while the models with axial ratio b/a in the range 
0.250 ^ b/a ^ 0.570 appeared stable. This paper extends that previous work to the study 
of a set of low angular momentum perfect elliptical disks. The minimal angular momentum 
cases allow us to study the ability of internal velocity dispersion to support an elliptic figure. 
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and forms a natural complement to the earlier work as the other bound of the whole class 
of perfect elhptic disks. As before, we tested for global stability by constructing a discrete, 
self-consistent model, loading it into an A^-body integrator and allowing it to evolve. 



2. Orbits in The Perfect Elliptic Disk 

Perfect elliptic disks are the two dimensional limiting case of the three dimensional 
perfect ellipsoids; their surface density S, derived from the perfect ellipsoid by integrating 
in z (see [Evans fc de Zeeuw 1992D , is given in Cartesian coordinates (x,?/) by 

(1 + m^y^i^ 0^ 

b < a, so that x is the major axis and y is the minor axis. The elliptic disk potential 
satisfies Stackel's criteria, implying that there are two integrals of motion [E and I2) and 
that the equations of motion are completely separable in confocal ellipsoidal coordinates 
( [Lynden-Bell 196^ ; |de Zeeuw 1985| ). The coordinates (A,/i) are the solutions to the 



quadratic equation 

2 2 

^— + ^ = 1, a<P<0, (2) 

where for the disk of equation (1) —a = a?, —j3 = b^. Curves of constant A are confocal 
ellipses aligned along the minor axis, with —a < A < 00, while curves of constant 
/i are hyperbolae, with — /3 < ^ < —a (see fig. |l]). All the curves share the foci 
X = 0, y = ±\//9 — a. When A ~ —a, the curves of constant A are highly elongated in the y 
direction; at large A, they become almost circular; A + a ~ r^. The curve /i = —(5 lies along 
the X-axis; as /i increases, the hyperbolae bend more sharply around the foci. 

The orbits in the elliptic disk potential divide into two families, box and loop orbits. 
The box orbits resemble Lissajous figures, combinations of independent oscillations in the 
X and y directions; they have no net angular momentum about the center of the potential. 
Loop orbits are ellipses or rosettes with a definite sense of rotation about the center of 
the potential. Because of the separability of the potential, the orbital bounding surfaces 
(extrema of A and ^) or turning points of the orbit are all lines of constant A or /i. Loop 
orbits are bounded by ellipses of constant A, given by the inner and outer turning points Ai 
and A2; a closed loop orbit is a curve of constant A. Box orbits are contained within an 
hyperbola of constant n and an ellipse of constant A, corresponding to the turning points 
/ii and A2. Only the loop orbits cross the minor axis outside the foci; within the foci, only 
the box orbits do so (see fig. |l|). An orbit may be characterized by its isolating integrals. 
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by its actions, or by its turning points. A more complete discussion, and transformations 
between the defining pairs can be found in LS, section 2. 



3. Constructing a self— consistent model 

Constructing a self-consistent model of a given potential amounts to finding a set 
of orbits specified by the values Ek, and l2k of E and I2 and a distribution weighting 
function Wk that together approximately reproduce the overall density distribution. The 
time-average densities Sorb along these individual orbits at any point (A, /i) must sum to 

^(•^'At) = H^orb(A,/i;-Efe,/2fe)wfc . (3) 

k 

There are two general approaches to the problem; we can either select a set of orbits by 
choosing the Ek and l2k, k = 1, ■ ■ ■ , N and then compute the distribution function weights 
Wk for each orbit, or we can choose the weights Wk, and then try to find a compatible set of 
orbits Ek and l2k- 

In the first method, the problem is substantially better constrained; if there are n 
orbits, already chosen, then we need only find n weights. In the second case, we have 
n weights, and are trying to find 2n orbit specifiers. The first approach has been used 
( [Schwarzschild 197^ ; Statler 1987; Teuben 1987; |de Zeeuw, Hunter fc Schwarzschild 1987. 



[hereafter ZHS| ) to demonstrate that solutions to the self-consistent problem do exist. If we 
wish to construct a discrete representation of a potential, for an A^-body simulation, then 
the second method has the advantage of allowing us to insist that each orbit contain an 
integral number of equal mass particles. LS showed that such problems could be solved 
successfully using the second approach; here we show that the first method can also be 
adapted for this purpose. Qualitative agreement between the two methods for the maximum 



angular momentum case has been shown in Levine (1995 



Along the minor axis, outside of the foci, the overall mass density is comprised solely 
of loop orbits; this provides a constraint that depends only on the selection of loop orbits. 
Loop orbits are also the only orbits with net angular momentum. This allows us to split 
the problem of choosing a set of orbits into two smaller problems. First, we select a set 
of loop orbits which add up to give the correct density in this region and also meet some 
additional criterion on the total angular momentum, and then we find a set of box orbits 
which contribute the rest of the mass needed. 
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3.1. Loop Orbit Selection 

If the outer turning point (A2) of a loop orbit is fixed, then as we move the inner 
turning point (Ai) outwards, the angular momentum of the orbit increases from a limit of 
zero for the marginal orbits (Ai = —a) up to a maximum which is reached at the closed 
loop orbit Ai = A2. Similarly, for a fixed value of Ai, the angular momentum of the orbit 
increases as we increase the value of A2. So, for maximum angular momentum, we choose 
the thinnest possible orbits (Ai = A2), and for the minimal angular momentum, the thickest 
possible (Ai <^ A2). 

In the maximum angular momentum case, the loop orbit population is easily and 
uniquely determined (ZHS; LS). The closed loops can never overlap each other, so at any 
point on the minor axis outside the foci, one and only one closed loop orbit contributes to 
the density at that point; finding the mass on each of the loop orbits is a one dimensional 
problem. In LS, we chose a set of n equally weighted orbits by integrating out from the 
minor axis; when the total integrated loop mass reached 1/n of the total loop orbit mass, 
we placed a single closed orbit on the mass weighted center ellipse of the elliptic annulus 
thus defined. The outer bound of the first orbit became the inner bound of the region 
represented by the next orbit, and the process was repeated until n orbits had been selected. 

Any other model must contain some thick loops, and very likely some of them will 
overlap each other, as well as the box orbits; finding a satisfactory set of loop orbits has 
become a two dimensional problem. We begin by dividing the minor axis into one 
dimensional cells, the inner-most cell boundary being the focal point (A = —a), and the 
outer-most being the outer hmit in A (Aout)- We create a library of N{N + l)/2 orbits 
whose inner and outer turning points lie on the boundaries of the cells (Acd)- For each outer 
turning point at Kdii = k) {k = 1, ■ ■ ■ , N) there are k inner turning points from \cd{i = 0) 
(close to the marginal orbit) out to Acd(^ = k — 1) (close to the closed loop orbit). From 
this library, we choose at most N orbits having non-zero weight, and compute weights Wcd 
for each member of the set of loop orbits. Because we must match the density in cells, 
and have A^(A^ -|- l)/2 possible orbits, the problem is under-constrained. The maximal and 
minimal angular momentum solutions are two cases where an angular momentum criterion 
coupled with the density matching restricts the possible solutions to a unique solution for a 
given library of orbits. 

To find a minimal angular momentum solution, we compute the weights sequentially, 
starting with the loops with their outer turning point in the outer-most cell. We take the 
loop closest to the marginal orbit (the "thickest" loop), and choose its weight so that it 
accounts for all of the mass in that cell. This fills the outermost cell with material from the 
orbit of lowest angular momentum which reachs that cell. The mass which must still be 
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placed in all the inner cells is then reduced by subtracting off the contribution of this orbit, 
while making sure that the density left in each cell is non-negative. If this is so, we move 
to the next outer-most cell, and repeat the process, and so step our way in towards the 
center. As we near the center, the contributions of the low angular momentum loops to the 
inner cells may fill some of them completely, before filling the outer-most loop. When this 
happens, we reduce the weight of the outermost orbit so that the remaining mass in each 
cell is non-negative, and then move the inner turning point out by one cell, and continue 
the process. We found that the loop orbit portion of the self-consistent model cannot be 
made up solely of marginal orbits; it is not possible to construct a model with exactly zero 
angular momentum. 

We have not proven formally that we have the minimum angular momentum case, but 
empirically we appear to be close to the minumum; the angular momenta of these models 
are ~ 10"^ times less than that of the maximum streaming models. As a check, we tried 
breaking up the solution space into an inner and an outer region at a A cell boundary, and 
solving first the outer and then the inner regions separately and adding up the solutions; no 
matter what boundary was chosen, the angular momentum was greater than in the original 
single region solution. In addition, the derivative of the angular momentum with respect to 
the inner turning point is quite steep, implying that we really want to push towards most 
marginal, as we have done. 

If the integrated linear density in the cells varies too greatly, we may not be able to 
find a positive solution. We separated the cell boundaries by uniform increments in their 
square roots, so that the cell divisions lie at 

AcdW = |v^=^ + (0^ - v^=^)^}' for i = (0, 1, • • • , AT). (4) 

(The small quantity S = 5 x 10~^ is added to the inner-most boundary to avoid the 
non-integrable singularity at A = —a.) The choice of square root spacing spreads the 
distribution function weights fairly evenly among the orbits. The problem can be solved 
with other spacings (we have used linear and exponential spacings as well though using an 
exponential cell spacing for the less centrally concentrated models did not always lead to a 
positive solution), but with much greater contrast in the orbit weights. 

As the axial ratio gets smaller, the gradient of the density along the minor axis (and 
overall) increases, making it desirable to use more cells to achieve a better approximation 
of the loop orbit density. If the number of cells on the minor axis is too small, then the 
discrete solution does not well approximate the continuous reality. Conversely, if the 
number is very large, then each orbit will end up being represented by a very small number 
of particles, and hence will itself not be well sampled. For b/a — 0.125, 0.570 and 0.910, 
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we computed solutions using between 100 and 1200 cells. We then chose the number of 
cells to be as small as possible and still give a good approximation to what appeared to be 
the limiting continuous solution. Figures ^ and ^ show respectively the relative weights of 
orbits as a function of turning points and the distribution of particles for 100 and 400 orbit 
library solutions. We decided to use 400 cells along the minor axis for all of the models; 
this provided a good compromise between good sampling of the linear density, and allowing 
high enough weights to permit us to sample each orbit reasonably well. 

3.2. Box Orbit Selection 

For the box orbits, the method of solution follows the method of ZHS; it is conceptually 
the same as for the loops, except that the grid is a full rectangle. First, the area of the disk 
is divided up into a grid in the coordinates (A,/i), with the grid points in both A and 
spaced by uniform increments in the square root (just as in the loop orbit case). For each 
grid cell, the coordinates of the corner with maximum A and /x are the turning points of the 
box orbits in the library. In each cell we compute the integrated surface density minus the 
sum of the surface densities already allocated to the loop orbits. For the box orbit with 
turning points in the outermost grid cell denoted (1, 1) in Fig. |^, the orbit's weight is chosen 
so as to supply all the mass in the cell. The contribution of this orbit is then subtracted 
from each of the inner cells, and the process repeated for the next outermost orbit and cell 
(denoted (1,2)), working inwards until weights have been computed for all the orbits. ZHS 
have previously shown that positive definite solutions do exist for the maximum streaming 
case; our experience is that solutions also exist for these minimum streaming models. 

We have chosen to use a library of orbits with 60 turning points in A and 30 in /x, 
making for 1800 orbits and cells. Figure ^ shows the difference in relative sampling between 
a library of 20 x 20 orbits and one with 60 x 30 orbits once they are populated with particles. 

4. Populating the Orbits 

We then place particles upon each orbit in a quiet manner, so as to minimize random 
noise in the initial conditions, and make it easier to watch for the growth of instabilities. 
The quiet distribution is our best approximation to a uniform distribution of particles 
throughout the phase space explored by an orbit. For a closed orbit, the solution is easily 
found (ISellwood 1983| ). For our space filling orbits in an integrable potential, we use a 
slightly modified version of the technique described in detail in LS. We place particles at 
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the intersection points of a grid on the torus in action-angle space which corresponds to 
each orbit. To find the starting positions and velocities of each particle, we construct a 
differential map from action-angle space to position and velocity coordinates. 

The weight of each orbit determines how many particles should be placed on the it. In 
all cases, we attempt to factor the integer nearest to this number into two factors ra, m as 
nearly equal as possible. If the nearest integer is less than 6, then we accept whatever pair 
of factors we compute. Otherwise, we also factor the second nearest integer, and choose the 
pair of factors that are most nearly equal; this helps to avoid problems when the nearest 
integer is prime, or has only two, very different, prime factors. We then use a grid with n 
lines evenly spaced in the angle 6\ and m lines in The overall difference between the 
desired and the computed number of particles is <^ 100 for models with 50,000 particles 
0.2% difference), and the model produces a good approximation of the overall potential. 



5. Model Setup and Integration 

All units from here on are expressed in terms where the gravitational constant G, the 
length scale a, and the total mass Mtotai of each model when integrated analytically out to 
infinity are all equal to 1. Because the perfect elliptic disks are formally infinite in extent, 
we have truncated the models, at r = 10a, for which all the models contain at least 90% of 
the total mass. Models were constructed for ellipticities ranging from b/a = 0.125 to 0.910 
following the prescription given above. These particle distributions were then loaded into 
an A^-body integrator and allowed to evolve under the infiuence of their own self-gravity. 
The fraction of the total mass within the truncation radius (Mtrunc); the axial ratio [b/a), 
the relative fraction of the model mass in the loop orbits, the angular momentum and the 
stability result for both the minimum and maximum cases of each of the models are given 
in Table |I[ 

We used a two dimensional polar-grid Fourier-transform A^-body code developed 
and kindly supplied by J. Sellwood (see ^ellwood 1981| , 1983 and [Miller 1976| for a more 



complete description). The grid was made up of 86 rings logarithmically spaced in radius 
with 100 grid points evenly spaced about them in azimuth. The grid was bounded by circles 
of radius r = 0.05a and r = 10.4a. The outermost cells had dimensions Ar = 0.64a by 
rA6 = 0.66a. Stars that crossed the outer boundary during the integration were discarded, 
while stars crossing into the central hole continued across it with constant velocity on a 



straight path, and were placed back on the grid at the next step ([Sellwood 1983|) . An 
explicit softening length of 0.05 was used for these integrations, and a small compensating 
radial force correction was added to insure that our initial models were close to virial 
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equilibrium. For an extended discussion of the details of the integration, and determination 
of the softening length see LS. 

The integration time step At was chosen to be less than 1/20**^ of the minimum 
period required for each of the angle variables to complete a circuit of 2tt radians on the 
action-angle torus for at least 90% of the orbits. We used At = 0.01 for all the models 
except b/a = 0.125 and 0.180, for which At = 0.002; the most elliptic models have higher 
central densities and velocities, requiring better time resolution. All of the models were 
run for at least 20 dimensionless time units, corresponding to 3-6 crossing times at the 
half-mass radius. This was enough time for gross instabilities to develop. For a number 
of the models we continued the simulations for 100 times units (15-30 half-mass crossing 
times) . 

At the beginning of the runs, the number of particles in the central hole was less than 
0.5%. In the most unstable models, several hundred particles fell off the grid within t = 20, 
while the stable models lost only tens of particles out of 50,000. The more elliptic models 
lost more particles, perhaps because of the higher radial velocities in their deeper central 
potentials. By t = 100, between 2% and 14% of the particles had left the grid. 



6. Results 

We constructed a set of minimal angular momentum models with ellipticities ranging 
from 0.125 to 0.910, highly elongated to almost circular. The minimal angular momentum 
models at both extremes of axial ratio appear to be unstable. For the nearly axisymmetric 
models (for b/a = 0.910 to 0.640, see fig. ^ and the dominant instability manifests itself 
as an increase in ellipticity, resembling the radial orbit instability seen in three dimensional 
simulations that have little rotational support (e.g. [Merritt fc Aguilar 1985| ; [Barnes, 



[Goodman fc Hut 1986|) . In the most elongated models (fig. |]) where b/a <^ 0.2, we see 
initially the beginnings of a bending instability, just as in the maximal streaming models, 
and similar to the bending seen in the prolate E9 model of Merritt & Hernquist (1991). 
This is followed later by a decrease in the overall ellipticity of the figure, and a growth 
in the power in the "lopsided" m = 1 mode similar to that seen in [Palmer fc Papaloizou 



(1990)1 and [Palmer, Papaloizou fc Allen (1990}] . Models with axial ratio ratio between 0.305 



and 0.570 appear to be stable over the duration of the simulations (fig. 0). 

Since the growth of an instability is not always apparent in the plots of particle 
position, we have also examined the behavior of the first six logarithmic spiral coefficients 
for these models. The growth of asymmetry in the m = 2 mode was defined in LS eq. [31] 
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as 

\A[m = 2,p = 0)\ 

this measures the spirahty which is inherent in growing modes ( |Lynden-Bell fc Ustriker 
|1967| ; LS, eq. [27] & [31]). Incipient spiral instabihties can show up here before they are 
clearly visible in the plots of particle positions. We call "stable" those models for which there 
is no change, above the noise level inherent in the particle discreteness, in the amplitude of 
the spiral harmonics. Figure |^ shows the change over the course of an A^-body integration 
of the dominant m = 2 log spiral mode, for the models of figures Bending-unstable 
models such as b/a = 0.125 grow asymmetrically at small p/m; this is most apparent as 
the s-shaped curve in the symmetry plot. The unstable nearly-round models which become 
more elliptic show a mostly symmetric growth in the m = 2 power. Stable models such 
ash/a = 0.305 don't budge. We have labeled the models with h/a = 0.250 and 0.640 as 
marginally stable: we see global instabihties of small power develop, but these very quickly 
saturate and die out. 



7. Discussion 

In this paper and in LS, we have constructed discrete self-consistent representations 
of the distribution functions of a range of perfect elliptic disks with minimal and maximal 
angular momentum. These models were then integrated forward in time using an iV-body 
integrator to see if they were stable. The nearly axisymmetric and the most elongated 
models were unstable. The perfect elliptic disks with moderate axial ratios appear to be 
stable in both the maximum and minimum streaming cases. 

In the maximum streaming case, the nearly axisymmetric models developed spiral and 
bar instabilities as expected, since their limiting case, a cold axisymmetric disk, is known 
to be violently unstable to spiral instabilities. In the minimal angular momentum case, 
nearly-round disks became more elliptical, in a manner very similar to the radial orbit 
instability of spherical systems. This is not too surprising, since the velocity distribution 
is anisotropic, with the radial velocity dispersion being substantially higher than the 
tangential dispersion, even in the very nearly axisymmetric models. This comes about 
because of the substantial presence of box and marginal loop orbits in the models. 

In both angular momentum extremes, the most elongated models developed a bending 
instability. The similarity in behavior is not very surprising given the decreasing importance 
of rotational support with increasing ellipticity in these models. P?remaine fc de Zeeuw| 
1987)1 have shown that the limiting case of the needle (6 0) is neutrally stable to bending. 
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while Merritt & Hernquist (1991) have demonstrated a bending instabihty in a very prolate 
(E9) system. It is thus not surprising that the most elliptic models should develop this 
instability. For the minimum angular momentum family, the instabilities change the shape 
of the disk towards a more moderate ellipticity. 

In the nearly axisymmetric disks, as the angular momentum is decreased from a 
maximum, we expect that the increasing velocity dispersion should help to stabilize against 
spiral instabilities. It appears likely that there is a stable region for nearly axisymmetric 
disks with values of Toomre's ( |1964| ) stability parameter Q which lies between the points 
Q ~ 2 and Q ~ 3 where the velocity dispersion has increased to the point of being 
able to support the disk against the spiral instability (fig |IU]). As the rotational support 
becomes negligible and the radial velocity dispersion increases, a radial-orbit instability 
develops; the disks with lower angular momentum become unstable to elliptical distortions 
when Tradiai /^tangential ^ 1-2 (after the discussion of [Fridman fc Polyachenko 1984| and 



Polyachenko 1987| for stability of spherical systems). We expect that there is a range of 
angular momentum between the two extremes for which the nearly round disks are stable. 
The moderately elliptical disks with maximum and minimum angular momentum appear to 
be stable, so we would anticipate that disks of similar ellipticity and intermediate angular 
momentum will also be stable. 

The stability of the two-dimensional models with moderate ellipticity gives us hope 
that the three dimensional perfect ellipsoids of intermediate triaxiality (which is probably 
the appropriate range for elliptical galaxies [Mihalas fc Binney 1981| ; |de Vaucouleurs,| 
de Vaucouleurs fc Cor win 1976| ; |de Zeeuw fc Franx 1991| ), will also be stable. It is known 



that some very flattened systems, such as the extreme oblate spheroids constructed from 
thin short-axis tube orbits ( [Merritt fc Stiavelli 1990| ), are unstable, but the simple fact 
that two longer axes are unequal is not likely to be the cause of further trouble. The three 
dimensional extension of this work will be interesting to see in light of the work of [Allen 



Palmer fc Papaloizou (1992)| showing that three dimensional systems with a small amount 



of rotational streaming are unstable to a tumbling bar instability, both when the models 
have largely radial orbits and when the orbits are mostly circular. 

The techniques developed in this work have laid the foundation for investigating the 
stability of three dimensional perfect ellipsoids, and indeed of any integrable potential. 
The methods for choosing orbits, whether simulated annealing (as in LS) or the marching 
scheme of ZHS, can be easily expanded to take account of a variety of possible cost terms 
related to angular momentum or line of sight velocities (e.g. [Kix 1997| ). The procedure of 
LS for generating a quiet start which minimizes random noise due to particle discreteness, 
can be carried over to any integrable system. This is potentially most useful in A^-body 
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studies which attempt to measure the growth rate of instabihties, because the detection of 
instabihties which are still in the linear regime is limited by particle noise. For example, 
Saha (1991)1 found that his linear stability theory was consistent with the results of A^-body 



simulations for highly unstable spherical systems, but predicted slow growing instabilities 



which could not be seen in the simulations because of particle noise. [Allen, Palmer 



Papaloizou (1990~)| have constructed an analytic potential-smoothing integration technique 
which decreases the \fN noise associated with binning and softening in iV-body codes, 
and permits better examination of the linear growth regime. Their method would also 
benefit from a quiet start, because the particle discreteness then makes a larger relative 
contribution to the noise. 

The authors would like to thank P. T. de Zeeuw for continued interest and for 
suggesting that we make sure that the ZHS and LS schemes agree, and J. Sellwood for 
graciously allowing us to use his A^-body integrator. This work has received support 
from grants 3739-E from the Consejo Nacional de Ciencia y Tecnologfa of Mexico, and 
NAGW-2769 from the National Aeronautics and Space Administration of the USA. 
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Table 1. Model Details 



Minimum Maximum 



h/a 


Mtrunc 


■Mtrunc 


10=^ X 


Stable^ 


Moop 
■Mtrunc 


10^ X 


Stable^ 


0.125 


0.935 


0.002 


0.012 


Bend 


0.01 


11.98 


Bend 


0.180 


0.934 


0.003 


0.025 


Bend 


0.02 


26.31 


Bend 


0.250 


0.932 


0.007 


0.051 


Marginal 


0.04 


50.84 


Stable 


0.305 


0.930 


0.010 


0.078 


Stable 


0.06 


76.68 


Stable 


0.370 


0.928 


0.015 


0.120 


Stable 


0.09 


114.80 


Stable 


0.440 


0.926 


0.022 


0.177 


Stable 


0.13 


165.27 


Stable 


0.470 


0.924 


0.026 


0.205 


Stable 


0.15 


189.75 


Stable 


0.570 


0.920 


0.040 


0.312 


Stable 


0.23 


284.66 


Stable 


0.640 


0.917 


0.054 


0.399 


Marginal 


0.30 


364.85 


Marginal 


0.715 


0.914 


0.074 


0.497 


Elliptic 


0.39 


463.27 


Spiral 


0.820 


0.909 


0.115 


0.629 


Elliptic 


0.55 


624.68 


Spiral 


0.910 


0.905 


0.176 


0.708 


Elliptic 


0.73 


784.92 


Spiral 



^Bend denotes bending into an 'S' shape, with decreasing b/a, Elliptic means unstable to 
increasing ellipticity, Spiral is unstable to forming a spiral. 



Loop Orbit 




Fig. 1. — Sample loop (top) and box (bottom) orbits showing the placement of particles upon 
each. The orbital bounding surfaces or turning points (ri and T2) are marked with dashed 
lines, and lines of constant angle have been plotted in the orbits to show the placement of 
each particle at an action-angle grid intersection. For the box orbit, the filled circles mark 
the intersections where particles are placed. For the loop orbit, for clarity, we have plotted 
only hnes of constant 9^; the intersections of these with lines of 9x are marked with filled 
triangles, and squares. The open squares mark the foci of the elliptic coordinates. 
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Fig. 2. — For a model with axis ratio b/a — 0.715, the top panel shows the turning points of 
the loop orbits with non-zero weight in the solution. The dashed line on the left is the line 
Ai = A2. The orbits with the outermost turning points are as close to the marginal orbit 
as the library contains. The middle panel shows the number of particles assigned to each 
orbit as a function of the outer turning point (A2). The 400 loop orbits contain a total of 
3760 particles. The bottom panel show the angular momentum on each orbit (solid line) 
and the cumulative angular momentum as a fraction of the total as a function of the orbit 
outer turning point (dashed line). The bulk of the angular momentum comes from the inner, 
non-marginal orbits. 
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:(a) b/a = 0.715 N^.u^^lOO 

[ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1(b) b/a = 0.715 N^,,i3 = 400 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1(c) b/a = 0.715 N^^jj^ = 20x20 

■ ■■■..J ■ .[.:■ 

~ . . ri.' 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1(d) b/a = 0.715 N^.^i,^ = 60x30 



Fig. 3. — For the model of figure 2, particle distributions on the loop orbits for libraries with 
(a) 100 and (b) 400 orbits, and on the box orbits for hbraries with N divisions in A and M 
divisions in where N x M are 20 x 20 (c) and 60 x 30 (d) orbits respectively. The small 
circles seen in (a) are a result of using too few orbits, with the result that some of the loop 
orbits have many more particles than the rest. The solution in (a) has 1814 particles, and 
that in (b) has 3760. 
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Fig. 4. — The ordering of the grid cells for the box orbit library, and their respective orbits 
in the elliptical space in which the equations of motion separate. Nine of the outermost 
cells and orbits are shown explicitly, with the orbit turning points labeled by filled squares, 
([a; = M, • • • , 1], 1) denote the marginal (thickest) box orbits. 
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b/a = 0.910 
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Fig. 5. — Positions {left column) and velocities (radial velocity: middle column; tangential 
velocity: right column) for 2500 of the 50,000 particles in the A'"-body integrations are shown 
at times 0, 5, 10, 15 and 20, for h/a = 0.910. 
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b/a = 0.640 
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Fig. 6. — Positions {left column) and velocities (radial velocity: middle column; tangential 
velocity: right column) for 2500 of the 50,000 particles in the A'"-body integrations are shown 
at times 0, 5, 10, 15 and 20, for h/a = 0.640. 



10 



(x, y) 



- 22 - 

b/a = 0.305 
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Fig. 7. — Positions {left column) and velocities (radial velocity: middle column; tangential 
velocity: right column) for 2500 of the 50,000 particles in the A'"-body integrations are shown 
at times 0, 5, 10, 15 and 20, for b/a = 0.305. 
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Fig. 8. — Positions {left column) and velocities (radial velocity: middle column; tangential 
velocity: right column) for 2500 of the 50,000 particles in the A'^-body integrations are shown 
at times 0, 5, 10, 15 and 20, for b/a = 0.125. 
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Fig. 9. — The left hand column shows the distribution of power with the tangent of the 
pitch angle in the m = 2 log spiral mode (In \A{m = 2,p)\) for models with b/a = 0.125, 
0.305, 0.640, and 0.910. The right hand column shows the growth of the asymmetry term 
A(m = 2,p) in the m = 2 mode. These are plotted at times {solid line), 10 {dotted line) 
and 20 {dashed line). 
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Fig. 10. — Disk angular momentum as a function of ellipticity for the maximum and 
minimal angular momentum disks. The solid line shows the angular momentum in the 
infinite maximum angular momentum disks, the squares the actual values in the truncated 
maximum angular momentum disk models, and the triangles the values for the minimal 
angular momentum disk models. Open and filled symbols represent unstable and stable 
models respectively. Values of Q for axisymmetric models are indicated along the right hand 
axis. Tr/Tt — 0.85 and 1.14 for Q = 2 and 3 respectively. 



